BuildNetwork

library(AnNet)
library(synaptome.db)
library(pander)
library(ggplot2)
library(ggrepel)

scale <- function(x, VALUE=NULL){

    x = as.numeric(as.vector(x))

    xmin <- min(x,na.rm=T)
    xmax <- max(x,na.rm=T)
   
    if( is.null(VALUE) ){

        x  <- x-xmin
        x  <- ifelse(!is.na(x), x/(xmax-xmin), NA) 

        return(x)
    }

    value = as.numeric(as.vector(value)[1])
    value = value-xmin
    value = ifelse(!is.na(value), value/(xmax-xmin), NA) 
    return(value)
}

Title AnNet: library for comprehensive network analysis

Authors: Colin Mclean, Anatoly Sorokin, J. Douglas Armstrong and Oksana Sorokina

#Introduction

Overview of capabilities

Build the network

The package allows building the network from the data frame, where the rows correspond to the edges of the graph, or for the list of nodes (genes), for which the information will be retrieved from the SynaptomeDB.

Build network for the node list extracted from SynaptomeDB


cid<-match('Presynaptic',getCompartments()$Name) # Let's get the ID for presynaptic compartment
cid
#> [1] 2
t<-getAllGenes4Compartment(cid) # Now we need to collect all the gene IDs for presinaptic  compartment
dim(t)
#> [1] 2304    8
head(t)
#> # A tibble: 6 × 8
#>   GeneID MGI       HumanEntrez MouseEntrez RatEntrez HumanName MouseName RatName
#>    <int> <chr>           <int>       <int>     <int> <chr>     <chr>     <chr>  
#> 1      1 MGI:1277…        1742       13385     29495 DLG4      Dlg4      Dlg4   
#> 2      2 MGI:88256         815       12322     25400 CAMK2A    Camk2a    Camk2a 
#> 3      3 MGI:96568        9118      226180     24503 INA       Ina       Ina    
#> 4      4 MGI:98388        6711       20742    305614 SPTBN1    Sptbn1    Sptbn1 
#> 5      5 MGI:88257         816       12323     24245 CAMK2B    Camk2b    Camk2b 
#> 6      6 MGI:1344…        1740       23859     64053 DLG2      Dlg2      Dlg2
gg<-buildFromSynaptomeByEntrez(t$HumanEntrez) # finally, build the graph using respecctive gene EntrezIDs as node IDs
summary(gg)
#> IGRAPH 42b7fb1 UN-- 1780 6620 -- 
#> + attr: name (v/c)

Annotate the nodes with node attributes

Gene name


gg<-annotateGeneNames(gg)
#> 'select()' returned 1:1 mapping between keys and columns
summary(gg)
#> IGRAPH 42b7fb1 UN-- 1780 6620 -- 
#> + attr: name (v/c), GeneName (v/c)
head(V(gg))
#> + 6/1780 vertices, named, from 42b7fb1:
#> [1] 273   6455  1759  1785  26052 2923
head(V(gg)$GeneName)
#> [1] "AMPH"   "SH3GL1" "DNM1"   "DNM2"   "DNM3"   "PDIA3"

Clustering

alg<- 'louvain'
gg<-calcClustering(gg,alg = alg)
{r cluster.graph, cache=TRUE}
gg<-calcAllClustering(gg)

Centrality

gg <- calcCentrality(gg)

Consensus matrix

Build consensus matrix for louvain clustering

conmat<-makeConsensusMatrix(gg,N=100,mask = 10,alg = alg,type = 2)
steps <- 100
Fn  <- ecdf(conmat[lower.tri(conmat)])
X<-seq(0,1,length.out=steps+1)
cdf<-Fn(X)
dt<-data.frame(cons=X,cdf=cdf)
ggplot(dt,aes(x=cons,y=cdf))+geom_line()+
      theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))

Cluster robustness

clrob<-getRobustness(gg,alg,conmat)
pander(clrob)
alg C Cn Crob CrobScaled
louvain 1 232 0.0999 1
louvain 2 204 0.09464 0.6492
louvain 3 295 0.09055 0.3769
louvain 4 194 0.09118 0.419
louvain 5 106 0.09491 0.6673
louvain 6 100 0.09529 0.6928
louvain 7 130 0.09105 0.4103
louvain 8 53 0.09707 0.8116
louvain 9 125 0.09104 0.4091
louvain 10 135 0.09146 0.4373
louvain 11 61 0.08832 0.2281
louvain 12 35 0.09178 0.4588
louvain 13 23 0.09693 0.8018
louvain 14 57 0.0849 0
louvain 15 30 0.09114 0.4161

Bridgeness

bm<-getBridgeness(gg,alg,conmat)
#> calculating Bridgeness for:  louvain
pander(head(bm))
ENTREZ.ID GENE.NAME BRIDGENESS.louvain
273 AMPH 0.0931096501093778
6455 SH3GL1 0.283225712276847
1759 DNM1 0.0932852071654777
1785 DNM2 0.45326719262814
26052 DNM3 0.0727508006268207
2923 PDIA3 0.581886549491273 
VIPs=c('8495','22999','8927','8573','26059','8497','27445','8499')
# VIPs=c('81876','10890','51552','5874','5862','11021','54734','5865','5864',
#            '9522','192683','10067','10396','9296','527','9114','537','535',
#            '528','51382','534','51606','523','80331','114569','127262','57084',
#            '57030','388662','6853','6854','8224','9900','9899','9145','9143',
#            '6855','132204','6857','127833','6861','529','526','140679','7781',
#            '81615','6844','6843')
indx   <- match(V(gg)$name,VIPs)
group <- ifelse( is.na(indx), 0,1)
MainDivSize <- 0.8
xmin        <- 0
xmax        <- 1
ymin        <- 0
ymax        <- 1
Xlab <- "Semilocal Centrality (SL)" 
Ylab <- "Bridgeness (B)"
X    <- as.numeric(igraph::get.vertex.attribute(gg,"SL",V(gg)))
X    <- scale(X)
Y       <- as.numeric(as.vector(bm[,dim(bm)[2]])) 
lbls <- ifelse(!is.na(indx),V(gg)$GeneName,"")
dt<-data.frame(X=X,Y=Y,vips=group,entres=V(gg)$name,name=V(gg)$GeneName)
dt_vips<-dt[dt$vips==1,]
dt_res<-dt[dt$vips==0,]
##--- baseColor of genes
baseColor="royalblue2"

##--- SPcolor, colour highlighting any 'specical' genes
SPColor="royalblue2"

##--- PSDColor, colour of core PSD95 interactor genes
PSDColor="magenta"

ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
    geom_point(data=dt_vips,
               aes(x=X,y=Y),colour=baseColor,size = 7,shape=15,show.legend=F)+
    geom_point(data=dt_res,
               aes(x=X,y=Y, alpha=(X*Y)), size = 3,shape=16,show.legend=F)+
    geom_label_repel(aes(label=as.vector(lbls)),
                     fontface='bold',color='black',fill='white',box.padding=0.1,
                     point.padding=NA,label.padding=0.15,segment.color='black',
                     force=1,size=rel(3.8),show.legend=F,max.overlaps=200)+
      labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
    scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) + 
    scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
    theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))+
        geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
    geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
#> Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

library(plotly)
#> 
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#> 
#>     last_plot
#> The following object is masked from 'package:latex2exp':
#> 
#>     TeX
#> The following object is masked from 'package:AnnotationDbi':
#> 
#>     select
#> The following object is masked from 'package:IRanges':
#> 
#>     slice
#> The following object is masked from 'package:S4Vectors':
#> 
#>     rename
#> The following object is masked from 'package:igraph':
#> 
#>     groups
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following object is masked from 'package:graphics':
#> 
#>     layout
p<-ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
    geom_point(data=dt_vips,
               aes(x=X,y=Y),colour=baseColor,shape=15,show.legend=F)+
    geom_point(data=dt_res,
               aes(x=X,y=Y, alpha=(X*Y)),shape=16,show.legend=F)+
      labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
    scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) + 
    scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
    theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
          panel.grid.minor = element_line(colour="grey40",size=0.1),
          panel.background = element_rect(fill="white"),
          panel.border = element_rect(linetype="solid",fill=NA))+
        geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
    geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
ggplotly(p)
#> Warning: `gather_()` was deprecated in tidyr 1.2.0.
#> Please use `gather()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.